{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Effective areas for IceCube Realtime Track Neutrino alerts\n",
    "\n",
    "Effective areas are a useful quantity for categorizing the overall sensitivity of a neutrino detector like IceCube, including triggering and event selection effects.  Given it's position at the South Pole, all variation in sensitivity can be quanitified in neutrino energy and declination.\n",
    "\n",
    "Additional information on calculation and use of effective areas can be found:\n",
    "\"All-sky Search for Time-integrated Neutrino Emission from Astrophysical Sources with 7 yr of IceCube Data\"\n",
    "Citation: M. G. Aartsen et al 2017 ApJ 835 151\n",
    "DOI: 10.3847/1538-4357/835/2/151\n",
    "\n",
    "We quantify the effective area separately for each of the IceCube internal event selections that contribute to the effective area. Details on each of these selection is in the supporting catalog paper:\n",
    "* Effa_gfugold.txt - GFU Gold seletion\n",
    "* Effa_ehegold.txt - EHE Gold selection\n",
    "* Effa_hesegold.txt - HESE Gold selection\n",
    "* Effa_gfubronze.txt - GFU Bronze selection\n",
    "* Effa_hesebronze.txt - HESE Bronze selection\n",
    "\n",
    "As the combination of these alerts with properly accounting for overlap of events in selections, we also provide the effective areas for the combined selection at the Gold level selections and the Gold + Bronze level selections, which serve as good representations of the overall sensitivty of the alert selections:\n",
    "* Effa_all_streams_gold_only.txt - Overall GOLD selection effective areas\n",
    "* Effa_all_streams_gold_bronze.txt - Overall GOLD+BRONZE selection effective areas\n",
    "\n",
    "These effective areas are binned in neutrino energy (10 TeV - 500 PeV) and in broad declination (stored internally in the files as local zenith angle) bins that reflect the major changes in effective area due to background rejection in the down-going direction and Earth absorption at higher energies for the up-going directions.\n",
    "* Declination boundaries (90 deg, 30 deg, 0 deg, -5 deg, -30 deg, -90 deg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot effective areas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "mpl.rc('font',family='serif',size=16)\n",
    "mpl.rc(\"text\",usetex = True)\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "#Generate Declination bins\n",
    "bins_coszen = [180, 120, 90, 85, 60, 0]\n",
    "coszen_min = bins_coszen[:-1]\n",
    "coszen_max = bins_coszen[1:]\n",
    "cmin = [a-90 for a in coszen_min]\n",
    "cmax = [a-90 for a in coszen_max]\n",
    "\n",
    "color_list = ['b', 'y', 'k', 'r', 'c']\n",
    "#Read effective area\n",
    "effA = np.genfromtxt('eff_area/Effa_all_streams_gold_bronze.txt')\n",
    "#effA = np.genfromtxt('eff_area/Effa_all_streams_gold_only.txt')\n",
    "\n",
    "\n",
    "ebins = effA[:,0]\n",
    "\n",
    "plt.figure(figsize = (8,8))\n",
    "ax1 = plt.subplot (1, 1, 1)\n",
    "\n",
    "for i,mx,mn in zip(range(1,6),cmax,cmin):\n",
    "    eA = effA[:,i]\n",
    "    ax1.plot(ebins,eA,drawstyle='steps-post', color=color_list[i-1], linewidth=2.5,\n",
    "                  label=r'${:.2f}^\\circ \\leq \\delta  < {:.2f}^\\circ$'.format((mx), (mn)))\n",
    "\n",
    "# Plot the overall stuff\n",
    "ax1.set_ylim(1e-1,5e4) \n",
    "ax1.set_xlim(1e4, 5e8)\n",
    "\n",
    "#plt.grid(True)\n",
    "ax1.set_yscale('log')\n",
    "ax1.set_xscale('log')\n",
    "ax1.set_ylabel(r'Effective Area (m$^{2}$)')\n",
    "ax1.set_xlabel('Neutrino Energy (GeV)')\n",
    "\n",
    "ax1.tick_params(axis='y', which='major', labelsize=16, length = 16, width = 1, direction = 'in', right = 'on', top = 'on')\n",
    "ax1.tick_params(axis='x', which='major', labelsize=16, length = 16, width = 1, direction = 'in', right = 'on', top = 'on', pad = 7)\n",
    "ax1.tick_params(axis='both', which='minor', labelsize=10, length = 8, width = 1, direction = 'in', right = 'on', top = 'on')\n",
    "ax1.legend(loc='upper left', ncol=1, fancybox=True)\n",
    "\n",
    "# Set grid lines\n",
    "for ymaj in ax1.yaxis.get_majorticklocs():\n",
    "        ax1.axhline (y=ymaj,ls=':', color='gray', alpha =0.75, linewidth=0.5)\n",
    "for xmaj in ax1.xaxis.get_majorticklocs():\n",
    "        ax1.axvline (x=xmaj,ls=':', color='gray', alpha =0.75, linewidth=0.5)\n",
    "\n",
    "plt.title('IceCube Track Alerts Gold + Bronze Effective Area')\n",
    "#plt.title('IceCube Track Alerts Gold Effective Area')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the effective area to estimate the number of events from the entire sky.\n",
    "\n",
    "To do this, the effective area over the 5 declination bands is flattened to a single all-sky declination band using\n",
    "the relative area of each band as a weight.\n",
    "\n",
    "Then a power-law astrophysical flux is convolved with the effective area to generate number of alerts expected per year"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "# Astrophysical muon neutrino flux found in:  https://iopscience.iop.org/article/10.3847/1538-4357/ac4d29\n",
    "index = -2.37\n",
    "flux =1.44e-18\n",
    "\n",
    "# Astrophysical muon neutrino flux found in: https://gcn.gsfc.nasa.gov/doc/IceCube_High_Energy_Neutrino_Track_Alerts_v2.pdf\n",
    "#index = -2.19\n",
    "#flux =0.9e-18\n",
    "\n",
    "\n",
    "days = 365.25\n",
    "\n",
    "\n",
    "Afile = np.genfromtxt('eff_area/Effa_all_streams_gold_bronze.txt')\n",
    "#Afile = np.genfromtxt('eff_area/Effa_all_streams_gold_only.txt')\n",
    "energy_bins = Afile[:,0]\n",
    "\n",
    "\n",
    "#flatten the A_eff from mutiple zenith bands\n",
    "bin_zenith_bound = [180, 120, 90, 85, 60, 0]\n",
    "bin_coszen_bound = np.cos(np.deg2rad(bin_zenith_bound))\n",
    "#print (bin_coszen_bound)\n",
    "A_eff_bins = Afile[:,1:]\n",
    "#print(A_eff_bins)\n",
    "all_zen_aeff = []\n",
    "for e_bin in range(0,len(energy_bins)):\n",
    "    this_bin_aeff = 0\n",
    "    norm=0\n",
    "    for item in range(1,len(bin_coszen_bound)):\n",
    "        this_bin_aeff += A_eff_bins[e_bin][item-1]*(bin_coszen_bound[item]-bin_coszen_bound[item-1])\n",
    "        norm += (bin_coszen_bound[item]-bin_coszen_bound[item-1])\n",
    "    all_zen_aeff.append(this_bin_aeff/norm)\n",
    "\n",
    "\n",
    "nus = 0 \n",
    "for i in range(1,len(energy_bins)): \n",
    "    #print (all_zen_aeff[i], energy_bins[i])\n",
    "    e_bin_width = energy_bins[i] - energy_bins[i-1] \n",
    "    nus += e_bin_width*all_zen_aeff[i] * 10000 * 86400 * days * 4*np.pi*flux *(energy_bins[i]/100000)**index\n",
    "print(nus)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
